Quantum Phase Transition in the Four-Spin Exchange Antiferromagnet 
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We study the S=l/2 Heisenberg antiferromagnet on a square lattice with nearest-neighbor and plaquette four- 
spin exchanges (introduced by A.W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).) This model undergoes a 
quantum phase transition from a spontaneously dimerized phase to Neel order at a critical coupling. We show 
that as the critical point is approached from the dimerized side, the system exhibits strong fluctuations in the 
dimer background, reflected in the presence of a low-energy singlet mode, with a simultaneous rise in the triplet 
quasiparticle density. We find that both singlet and triplet modes of high density condense at the transition, 
signaling restoration of lattice symmetry. In our approach, which goes beyond mean-field theory in terms of the 
triplet excitations, the transition appears sharp; however since our method breaks down near the critical point, 
we argue that we cannot make a definite conclusion regarding the order of the transition. 

PACS numbers: 75.10.Jm, 75.3().Kz, 75.50.Ee 



I. INTRODUCTION 

Problems related to quantum criticality in quantum spin 
systems are of both fundamental and practical importancai. 
Numerous materials, such as Mott insulators, exhibit either 
antiferromagnetic (Neel) order or quantum disordered (spin 
gapped) ground state depending on the distribution of Heisen- 
berg exchange couplings and geometry. External perturba- 
tions (such as doping or frustration) can also cause quantum 
transitions between these phases. Systems with spin 1/2 are 
indeed the most interesting as they are the most susceptible to 
such transitions. It is well understood that the quantum transi- 
tion between a quantum disordered and a Neel phase is in the 
0(3) universality classi, where a triplet state condenses at the 
quantum critical point (QCP). 

A recent exciting development in our theoretical under- 
standing of QCPs originated from the proposal that if the 
quantum disordered (QD) phase spontaneously breaks lattice 
symmetries {e.g. is characterized by spontaneous dimer or- 
der), and the transition is of second order, then exactly at 
the QCP spinon deconfinement occurs, i.e. the excitations 
are fractionalized^. It is assumed that the Hamiltonian itself 
does not break the lattice symmetries {i.e. does not have "triv- 
ial" dimer order caused by some exchanges being stronger 
than the others). We use the terms "dimer order" and "va- 
lence bond solid (VBS) order" interchangeably. It is expected 
that the dimer order vanishes exactly at the point where Neel 
order appears, i.e. there is no coexistence between the two 
phases. Deconfinement thus is intimately related to disap- 
pearance of VBS order; indeed if the latter persisted in the 
Neel phase it would be impossible to isolate a spinon, as 
"pairing" would always take place. Spontaneous VBS or- 
der driven by frustration has been a common theme in quan- 
tum antiferromagnetism^., although its presence and the na- 
ture of criticality in specific models, such as the 2D square- 
lattice frustrated Heisenberg antifeiTomagnet, is still some- 
what controversial^. It would be particularly useful to apply 
unbiased numerical approaches, such as the Quantum Monte 
Carlo (QMC) method, to study frustrated spin models; how- 
ever due to the fatal "sign" problem^, frustrated Heisenberg 
systems are beyond the QMC reach. 



In a recent study, the QMC method was applied to a four- 
spin exchange quantum spin model without frustration, which 
was shown to exhibit columnar dimer VBS order and a mag- 
netically ordered phase with a deconfined QCP separating 
them^. These conclusions were later confirmed by further 
QMC studies^. Extensions of the model, which include for 
example additional (six-spin) interactions, provide additional 
support for a continuous QCP**. A different VBS pattern (pla- 
quette order) was also proposed for the four-spin exchange 
model^. At the same time, the nature of the quantum phase 
transition was challenged in Refs. where arguments 

were given that the transition is in fact of (weakly) first order. 

It is the objective of the present work to study the Sand- 
vik model^, by approaching the quantum transition from the 
dimer VBS phase. Our approach uses as a starting point a 
symmetry broken state {i.e. one out of four degenerate VBS 
configurations), and we thus must search for signatures that 
the system attempts to restore the lattice symmetry at the QCP. 
Even though full restoration is impossible within the present 
framework, we find a QCP characterized by condensation of 
triplet modes of high density; this is in contrast to the conven- 
tional situation when the condensing particles are in the dilute 
Bose gas limit. The high density itself is due to the presence 
of a singlet mode that condenses at the QCP, and reflects the 
strong fluctuations of the background dimer order The above 
effects lead to the vanishing of the VBS order parameter; at 
the same time our method, which accounts for the strong fluc- 
tuations, leads to a rather sharp phase transition. It appears 
that we cannot draw a definite conclusion about the order of 
the transition because in the vicinity of the QCP the triplon 
density increases uncontrollably, suggesting that other states 
(such a plaquette states and larger clusters) are strongly ad- 
mixed into the ground state. This is generally expected in a 
situation where the lattice symmetry is restored at the quan- 
tum critical point. 

The model under consideration is 

H = jY^ Sa.Sfa -K Y^i^a- Sb){S, ■ Sd), (1) 

(a.b) a,b,c,d 

where J > 0, i^T > 0, and all spins are S* = 1/2. Consider 
the numbers 1,2,3,4 in Fig. [T] The summation in the four- 
spin term is over indexes (a, b) — (1, 2), (c, d) = (3, 4) and 
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FIG. 1: (Color online) Dimer pattern in the quantum disordered 
(VBS) phase, KjJ > {K/J)c. 



{a, b) = (1, 4), (c, d) — (2, 3) on a given plaquette, and then 
summation is made over all plaquettes'^. The range of pa- 
rameters explored in Ref. is K/ J < 2, and the QCP is at 
{K/ J)c ~ 1.85. Our coupling notation is slightly different 
from the one used in Refs.|6',7]; the coupling K is related to 
the parameter Q*"-^. via K = Q/{1 + Q/{2J)), and the criti- 
cal point in that notation is {Q/J)c ~ 25. The dimerization 
pattern is proposed to be of the "columnar" type, as shown in 
Fig-in Four such configurations exist. We will assume a con- 
figuration of this type, will show that it is stable at K/J 3> 1, 
and will then search for an instability towards the Neel state 
as K/J decreases. 

The rest of the paper is organized as followed. In Sec- 
tion II we present results based on the mean-field approach in 
terms of the dimer (triplon) operators. In Section III we extend 
our treatment beyond mean-field, and even further in Section 
IV, where we also take into account low-energy singlet two- 
triplon excitations. Section V contains our conclusions. 



II. MEAN-FIELD TREATMENT 

We start by rewriting Eq. ([T]l in the the bond-operator 
representation'^, where on a dimer i, the two spins forming it 
are expressed as: — ^{^sltia±tl^Si~i£ai3ftlpti^), and 
Sj , ij^, a — x,y, z create a singlet and triplet of states. We re- 
fer to the triplet (S=l) quasiparticle, tt^, as "triplon". The 
bold indexes i, j, m, 1 label the dimers (see Fig. [T]i. Summa- 
tion over repeated Greek indexes is assumed, unless indicated 
otherwise. 

The hard-core constraint, Sj Si + tl^Ua = 1, must be en- 
forced on every site, which at the mean-field (MF) level can be 
done by introducing a term in the Hamiltonian, — /i + 
tl^tia — 1). Then fi and the (condensed) singlet amplitude 
s = (si), are determined by the MF equations'^. We obtain at 
the quadratic level, in momentum representation: 

H2=J2 l^ktL^k. + ^ (C^lk. + h.c.) I (2) 



where 

= + s4E(k) , 

= -{J/2)cosk^ + {J±K/4)cosky . (3) 

The four-spin interaction from ([T]) acting between two 
dimers (e.g. i, j in Fig.[T]i contributes to the "on-site" gap and 
hopping (^jjT) via Ak, as well as to the quantum fluctuations 
term The part involving four dimers has been split in a 
mean-field fashion, leading to the Hartree-Fock self-energy 

— T,{h)/K — 2'Tix cos -\-2Yiy cos ky + Yi^y cos k^ cos ky , 

(4) 

with 

a 

where i, m are neighboring dimers in the x (horizontal) direc- 
tion (Fig.[ril, and similarly for the y and the diagonal contribu- 
tions. The triplon dispersion is uj (k) = V^k - B'i, and has 
a minimum at the Neel ordering wave-vector \^af = (0, tt) 
(since we work on a dimerized lattice). The ground state en- 
ergy is then easily computed, 

Egs ^Eo + {H2) , (6) 

where 

Eq/N = ^^{Js^+Ks^)+fi{~s^ + l)+ (7) 

3Ks^(^j:l + j:l + ^^%y (8) 

and 

k 

The mean-field equations require a numerical minimization 
of Egs with respect to the parameters {/i, s, E^;, Sy, S^^y}. 
This amounts to the self-consistent Hartree-Fock approxima- 
tion for I](k). The result for the triplon gap A — uiCkAp) is 
presented in Fig. |2] (black curve). 

The MF result {K/J)c ~ 0.6 substantially underestimates 
the location of the critical point, compared to the the QMC 
calculations, where {K/J)c ~ 1.85^'^. Interestingly, if one 
solves the MF equations ignoring both the hard core and the 
E(k), one finds (A7 J)c = 1. Physically, in the full MF, the 
hard core contribution increases the gap (and hence the stabil- 
ity of the dimer phase) while at the same time suppressing the 
antiferromagnetic fluctuations (which favor the Neel state). 

We also note that a recent (hierarchical) MF treatment 
based on the plaquette ground state also underestimates very 
strongly the QCP location ({K/J)c ~ 1^), similarly to our 
result. In our view this means that both mean field approaches 
are not sufficient to attack the present problem, where fluctu- 
ations are apparently very strong. We choose to accept that 
the numerical QMC result gives the most accurate determina- 
tion of the QCP, and therefore in what follows we extend our 
treatment in several directions beyond mean-field theory. 
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FIG. 2: (Color online) Triplon excitation gap A — uj(kAF) in vari- 
ous approximations. The point A ^ corresponds to transition to 
the Neel phase. 



III. BEYOND MEAN-FIELD: THE DILUTE TRIPLON GAS 
APPROXIMATION 

A more accurate treatment of fluctuations is possible by 
taking into account the hard-core constraint beyond mean- 
field. One can set the singlet amplitude s = 1 in the pre- 
vious formulas, but introduce an infinite on-site repulsion be- 
tween the triplons, U J^i ap ^Ii^/3i^/3i^aij ^ oo. As long as 
the triplon density (determined by the quantum fluctuations) is 
low, an infinite repulsion corresponds to a finite scattering am- 
plitude between excitations and can be calculated by resum- 
ming ladder diagrams for the scattering vertex'"*. This leads 
to the effective triplon-triplon vertex r(k, w) which was pre- 
viously calculatedi^: 

r-^(k,c.)^V ^ ^ "^"^-'^ ^ — + 1 j. 

(10) 

This vertex in turn affects the triplon dispersion via (what we 
call) the Brueckner self-energyi^: 

EB(k,c^) =4^t.2r(k + q,c^-t^(q)). (11) 
q 

The corresponding parameters in the quadratic Hamiltonian 
(O in this case are 

= J + 2if(l-4nt/3)+ek +S(k) + SB(k,0), 
= ^i[ + S(k). (12) 

The Bogolubov coefficients are defined in the usual way 
ul. = 1/2 + Ak/(2w(k)) = 1 + Wk- The various terms 
in I](k) can be expressed through them; for example Y^x = 
X]k(^k + ^^k^k) cosfc^:, and so on. The density of triplons 
is rit = {t\^tia) = SX^k'^k- In addition, the renormaliza- 
tion of the quasiparticle residue, Z^^ = 1 — (9S]B(k, 0)/duj, 



FIG. 3: Renormalization of quantum fluctuations by resummation of 
a ladder series, with l ll3b at the vertices. 

implies the replacement Uk — > \/^itk, vy^ VZ^vy^ in 
all the formulas'^, and the renormalized spectrum cj(k) = 

An iterative numerical evaluation of the spectrum using 
the above equations, which amounts to solution of the Dyson 
equation, leads to the result shown in Fig.|2](blue curve). The 
above approach appears to be well justified since the quasi- 
particle density rif < 0.1. The resulting critical point is still in 
the "weak-coupling" regime K/J < 1, with about 100% de- 
viation from the QMC result {{K/ J)c ~ 1.85). This suggests 
that the on-site triplon fluctuations are not the dominant cause 
for the disagreement with the QMC results; thus we proceed 
to include two-particle fluctuations (in the triplon language), 
which amounts to including dimer-dimer correlations. 

IV. STRONG FLUCTUATIONS IN THE SINGLET 
BACKGROUND: QCP BEYOND THE DILUTE TRIPLON 
GAS APPROXIMATION 

It is clear that "non-perturbative" effects are responsible for 
driving the QCP towards the "strong-coupling" region K/ J ^ 
2. To proceed we make two improvements to the previous 
low-density, weak-coupling theory. 

First, we take into account fluctuations in the singlet back- 
ground, i.e. the manifold on which the triplons are built and 
interact. The main effect originates from the action of the 
four-spin X-term from ([T]i on two dimers, e.g. i, j in Fig. [T] 
Part of this action has led to the on-site gap 2K in ( fT2l l. fa- 
voring dimerization. However, a strong attraction between the 
two dimers is also present, since the K-Xem\ is symmetric with 
respect to the index pair exchange (1,2)(3,4) ^ (1,4)(2,3), 
leading to a "plaquettization" tendency as well. In the triplon 
language this is manifested by formation of bound states of 
two triplons, due to their nearest-neighbor interactions 

Hi.y = {7itli4j*/3i*"j + 72tliilji/3ii/3j 

+ l^tlAs^o^itpi] , (13) 
K J K J 5K 

We also checked that on the perturbative (Hartree-Fock) level, 
the effect of this term on equations (O and (fT2] i was negligible 
(and we did not write it explicitly). 

An intuitive way of taking into account the effect of two- 
triplon bound states (with total spin S=0) on the one-triplon 
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FIG. 4: (Color online) (a.) Singlet bound state energy (black), 
binding energy e = 2A — Es (blue), and the triplon gap A (red), 
(b.) Triplon density nt. (c.) Dimer order parameters. Dashed parts 
of the lines represent points corresponding to rapid growth of the 
quasiparticle density. 



spectrum, is to work in the "local" approximation. This means 
effectively neglecting the triplon dispersion and directly eval- 
uating the ladder series that renormalizes the quantum fluctu- 
ation term Bk in (|2]l, corresponding to emission of a pair of 
triplons with zero total momentum. This is illustrated graphi- 
cally in Fig. |3] with the result 



Bl- = — — cosfc 




cos ky + S(k) 



7 = 7i + 372 + 73 = - J - -K, 



(14) 



where 7 is the effective attraction of two triplons with total 
S' = 0, and 



is the energy of two (non-interacting) triplons on adjacent 
sites. This calculation is justified for K/ J 3> 1 and leads 
to an increase of the quantum fluctuations, and from there to 
almost doubling of the triplon density nt (see Fig. |4] below). 
It contributes significantly to the shift of the QCP. 

We can go beyond the "local" approximation by solving 
the Bethe-Salpeter equation for the bound state, formed due 
to the attraction (fT3T l. and taking into account the full triplon 
dispersion. The equation for the singlet bound state energy 
Es{Q), corresponding to total pair momentum Q is 



27 



cj(Q/2 + q)-tj(Q/2-q) 



(16) 



AE = 2J+—K 
4 



(15) 



Here we have, for simplicity, written only the main contribu- 
tion to pairing (Eq. (fT3b ) in the limit K/ J ^ 1, and have 
neglected the on-site repulsion (which leads to slightly di- 
minished pairing), as well as small pairing due to the ex- 
change J from dimers in the a;-direction on Fig. [T] It is 
easily seen that the lowest energy corresponds to Q = 0; 
we define from now on Eg = Es{Q — 0). The bind- 
ing energy is e = 2 A — Eg, where A is the one-particle 
gap. The bound state wave-function corresponding to Eg is 
l*> - Ec,ij,,„*?.e'9''^'"-''4<ij|0). In the "local" limit 
(nearest-neighbor pairing), ^fq^^ = y/2 cos qy. 

Second, we have made subtle changes to the resumma- 
tion procedure concerning the quasiparticle renormalization 
Z, based on both formal and physical grounds. On the one 
hand it is clear that in the Brueckner approximation (Eq. (fTTTl). 
where the self energy is linear in the density (E^ oc nt), the 
dependence of the vertex F^^ on density is beyond the accu- 
racy of the calculation, meaning one can put Uq = l,Vq = 
in ( [Tol l, instead of determining them self-consistently. This 
leads to a decreased influence of the hard-core (which fa- 
vors the dimer state) on the Hartree-Fock self-energy E(k) 
from (|4|i (which favors the Neel state). It is indeed the mutual 
interplay between E^ > and E(k) < 0, that determines the 
exact location of the QCP in the course of the Dyson's equa- 
tion iterative solution. While in the "weak-coupling" regime 
K/ J < 1, Es always dominates, in the "strong-coupling" re- 
gion K/J > 2, E(k) starts playing a significant role, since 
parametrically E cx Knt- It is physically consistent that in 
the region where singlet fluctuations in the dimer background 
are strong, the hard-core effect is less important, i.e. in ef- 
fect the kinematic hard-core constraint is "relaxed". We also 
observe that in typical models with QCP driven by explicit 
dimerization, such as the bilayer model, the described differ- 
ence in approximation schemes makes a very small difference 
on the location of the QCP'^, since those models are always 
in the "weak-coupling" regime, dominated by the hard-core 
repulsion of excitations on a non-fluctuating dimer configura- 
tion. The purpose of the above rather technical diversion is to 
emphasize that care has been taken to take into account as ac- 
curately as possible the effect of the (low-energy) two-particle 
spectrum on the one-particle triplon gap. 

Our results are summarized in Fig.|4]and Fig.|2](red line) for 
the gap. The critical point is shifted towards (K/ J)c « 2.16 
(in much better agreement with QMC data), with a very 
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strong increase of the density towards Kc- This translates 
into a decrease of the dimer order, as measured by the two 
dimer order parameters that we compute from the expres- 
sions: - |(S3 • S4) - (S5 • 84)1 = I - I + nt + 
Dy = |(S3-S4)-(Si-S4)| = |-| + nt-fSy|. The spins are 
labeled as in Fig.[T] The singlet bound state energy (0) also 
tends towards zero at the QCP, with the corresponding binding 
energy remaining quite large e/ J w 1. All these effects point 
towards a tendency of the system to restore the lattice sym- 
metry, although it is certainly clear that as the critical point is 
approached, our approximation scheme (low density of quasi- 
particles) breaks down (dashed lines on figures). We should 
point out that the sharpness of variation near Kc is not due to 
divergence in any of the self-energies but is a result of rapid 
cancellation at high orders {i.e. iterations in the Dyson equa- 
tion). In fact cutting off our iterative procedure at finite order 
gives a smooth curve, suggesting that additional classes of di- 
agrams become important (although in practice their classifi- 
cation is an insurmountable task). The merger of singlet and 
triplet modes, which we find near the QCP, in principle reflects 
a tendency towards quasiparticle fractionalization (spinon de- 
confinement) and is also found in the ID Heisenberg chain 
with frustration'^, where spinons are always deconfined. 

Since we are now dealing with a situation where the density 
is not very small nt « 0.2, it is prudent to check how the next 
order in the density may affect the above results. For example 
at second order in the density, the self-energy S(k) changes 
by amount 5E(k), i.e. one has to add this contribution to the 
right hand side of Eq.(fT2]). namely: 

Ak^Ak + 5S(k), Sk^Bk + (5S(k). (17) 

We have found 

5E(k) = -2K{Ql- P'^)coBk^ + 
+2K{Ql~ P^) cos ky - 
-K{Qly- P^y) cos k^ cos ky, (18) 

and the following definitions are used: 

a k 

Qx = (1/3) ^(^iaW) = ^MkWkCOsfc^r, (19) 
a k 

and similarly for the other directions, for example Py = 
Sk ^k '^'^^ — J2k ^k After includ- 

ing these expressions in our numerical iterative procedure, we 
have found that the QCP is shifted by a very small amount, 
and the overall picture, as summarized in Fig. |4] and Fig. |2] 
(red line), still stands. 



V. CONCLUSIONS 

In conclusion, we have shown that the QCP between the 
Neel and the dimer state in the model ([T]i is of unconven- 



tional nature, in the sense that it is characterized by the pres- 
ence of both triplet and singlet low-energy modes. Near the 
QCP, whose location {{K/ J)c ~ 2.16) we find in fairly good 
agreement with recent QMC studies, the system exhibits: (1.) 
Strong rise of the triplon excitation density, due to increased 
quantum fluctuations, (2.) Corresponding strong decrease 
(and ultimately vanishing) of the dimer order at the QCP (3.) 
Vanishing of a singlet energy scale, related to the destruction 
of the dimer "columns" in Fig.[T] The above effects are all re- 
lated and influence strongly one another, ultimately meaning 
that the QCP reflects strong fluctuations and can not be de- 
scribed in a mean-field theory framework. These results also 
suggest a desire of the system to restore the lattice symmetry 
at the QCP, as found in the QMC studies*'. 

At the same time all our improvements beyond mean-field 
theory have also resulted in a very sharp transition, which ap- 
pears to be first order. However in our view our approach 
is not capable of addressing correctly the issue of the order 
of the phase transition, basically because once we take the 
strong (inter) dimer fluctuations in to account, the triplon den- 
sity starts rising quickly beyond control. This is in a certain 
sense natural in a situation where the system wants to restore 
the lattice symmetry at the QCP and thus the ground state 
acquires strong admixture of plaquette, etc. fluctuations as 
the dimers begin to "disappear." This is also manifested in the 
fact that our procedure is sensitive to the number of iterations 
in the Dyson equation; all presented results are for an "infi- 
nite" number of iterations, so that a fixed point is reached, 
but cutting off the procedure results in a smoother behavior 
and a shift of the QCP, which becomes iteration dependent. 
We have not previously encountered such volatile behavior in 
any other spin model with a dimer to magnetic order transi- 
tion. Since iterations translate into accounting of more and 
more fluctuations, the sensitivity of the results seems to mean 
that the situation starts spiraling out of control near the QCP, 
quite likely because classes of fluctuations become important 
that are not included in the dimer description, such as longer 
range correlations, etc. All this suggests that the triplon quasi- 
particle description breaks down near the QCP which indeed 
appears natural in a model where spinon deconfinement is ex- 
pected to take place at the QCP^. On the other hand, if we 
put aside the arguments that our approach is not reliable near 
the QCP, the natural conclusion would be that the transition is 
first order. 
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